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We use a mesoscopic computer simulation method to study the interplay between hydrodynamic 
and Brownian fluctuations during steady-state sedimentation of hard sphere particles for Peclet 
numbers (Pe) ranging from 0.1 to 15. Even when the hydrodynamic interactions are an order 
of magnitude weaker than Brownian forces, they still induce backflow effects that dominate the 
reduction of the average sedimentation velocity with increasing particle packing fraction. Velocity 
fluctuations, on the other hand, begin to show nonequilibrium hydrodynamic character for Pe I 1. 
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Exactly what happens when a collection of parti- 
cles sediments through a viscous solvent is a very sim- 
ple question to pose, but a remarkably difficult one to 
answer^, Q- In a classic tour de force, Batchelor^ 
showed that the average sedimentation velocity v s of hard 
spheres (HS) of hydrodynamic radius a has a lowest or- 
der correction v s = v°(l — 6.55</>) where is the HS 
volume fraction and v° s is the Stokes sedimentation ve- 
locity of a single sphere 0. This substantial correction to 
v® with volume fraction is caused by many-body hydro- 
dynamic backflow effects that greatly complicate efforts 
to extend the Bachelor result to higher order in </>. Even 
less is understood about the velocity fluctuations around 
the average, 5v = v — v s . Using straightforward physical 
arguments, Caffish and Luke|4| predicted that for sedi- 
mentation these should diverge as ((<5v) 2 ) ~ L, where L 
is the smallest container size. This surprising result pro- 
voked a flurry of experimental and theoretical studies (see 
ref.Q for a review). Although it is agreed that hydro- 
dynamic velocity fluctuations are relatively large, there 
is no consensus on the reasons (if any) for the purported 
breakdown of the Cafiish-Luke argument at large L. 

In addition to the fundamental interest and myriad ap- 
plications of sedimentation itself, researchers have been 
motivated to investigate this problem because of its rele- 
vance to non-equilibrium statistical mechanics. Recent 
studies in this line include a theoretical prediction of 
a continuous nonequilibrium noise-driven phase transi- 
tion between screened and unscreened phases |5j, and an 
experimental study predicting a noise-induced effective 
"temperature" that could aid in developing an ensemble 
based statistical mechanics for driven systems 6] . 

Most theoretical studies of sedimentation have focused 
on the limit where Brownian forces are negligible, and 
only hydrodynamic interactions (HI) contribute. In other 
words, the dimcnsionless Peclet number 



which measures the relative strength of HI and Brown- 
ian forces, was assumed to be infinite. Here D co i is the 



equilibrium self-diffusion constant of the particles. When 
the gravitational energy gained by a particle sedimenting 
over a distance of one radius a is equal to the reduced 
temperature fc^T, then Pe = 10, a criterion used to de- 
fine the start of the colloidal regime^]. Sedimentation at 
Pe < 1 has many important applications for colloidal dis- 
persions, as well as for centrifugal diagnostic techniques 
commonly used for biological macromoleculesQ. 

In this Letter, we employ a recently proposed meso- 
scopic simulation method |?J to investigate steady state 
sedimentation at finite Pe, where Brownian and HI both 
contribute to velocity fluctuations. To our knowledge, 
this problem has not been investigated in detail before. 
For all Pe studied, we find that the average sedimentation 
velocity is completely dominated by HI, even when they 
are much smaller than the Brownian forces. On the other 
hand, we argue that short time velocity fluctuations are 
dominated by Brownian forces up to surprisingly large 
Pe, while long time fluctuations have predominantly hy- 
drodynamic character even at moderate Pe. 

To perform the simulations, we adapt Stochastic Rota- 
tion Dynamics (SRD) [t| to the problem of sedimenting 
HS. SRD is a particle based method similar in spirit to 
the lattice Boltzmann model (LBV which has been ex- 
tensively applied to sedimentation^. In contrast to LB, 
it naturally includes Brownian noise (see however |9(). 
In SRD a fluid is represented by Nf ideal particles of 
mass rrif. After propagating the particles for a time At c , 
the system is partitioned into cubic cells of volume giq. 
The velocities relative to the center of mass velocity of 
each separate cell are rotated over a fixed angle around 
a random axis. This procedure conserves mass, momen- 
tum, and energy, and yields the correct hydrodynamic 
(Navicr Stokes) equations, including the effect of thermal 
noise Q • The fluid particles only interact with each other 
through the rotation procedure, which can be viewed as a 
coarse-graining of particle collisions over time and space. 
For this reason, the particles should not be interpreted as 
individual molecules but rather as a Navier Stokes solver 
that naturally includes Brownian noise. 

The colloid-colloid (cc) and colloid-fluid (cf) interac- 
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tions are modeled by a repulsive potential: f3V C i(r) = 
10 [{a cl /r) 2n - {a cl /r) n + 1/4] (r < 2 1 ' n <j cl ). For 
Kjc( r )> n — 24, which should approximate HS behav- 
ior. For V c f(r) we set n = 6 and a c f = 0.465cr cc , slightly 
below half the colloid diameter c cc , which allows for lu- 
brication. These potentials result in a hydrodynamic ra- 
dius a ss 0.8er c /. Colloid-colloid and colloid- fluid forces 
are integrated with a standard velocity Verlet molecular 
dynamics integrator with a time step t = \ A^|Toj. 

We now briefly discuss our choice of SRD parameters; 
a more detailed account will be published elsewhere |Tl| . 
The kinematic viscosity v = r]f/pf, where 77/ is the vis- 
cosity and pf the mass density of the fluid, is an impor- 
tant parameter because it sets the timescale over which 
the momentum (vorticity) diffuses away. In dimension- 
less form it is desirable to have v/D co i > Sc = v/Df 3> 1, 
where Sc is the Schmidt number and D t the self-diffusion 
constant of the fluid particles. Since D co i < Df, the 
first inequality is always satisfied. When Sc ~ 1 mo- 
mentum diffusion is dominated by mass diffusion, as in 
a gas. If Sc » 1 the fluid is liquid like. Since the 
SRD particles could be viewed as collections of individ- 
ual molecules, the Schmidt number of an SRD fluid is 
smaller than the real fluid it represents. Nevertheless, to 
model a liquid-like system, we choose a relatively small 
collision interval (At c = 0.1ao(m//fcsT) 1 / 2 ), leading to 
Sc Ri 5. To prevent compressibility effects, the gravita- 
tional field g was limited such that 0.0067 < v s /cf < 0.1, 
where cj — ^/2fcsT '/m/ is the speed of sound in the 
fluid. Finally, to avoid large inertial effects, the particle 
Reynolds number Re = v s ajv should be <C 1, as in real 
suspensions^. Inevitably there will be a compromise be- 
tween computational efficiency and low Re. In our work 
0.0016 < Re < 0.24, depending on Pe, which is similar 
to the choice made for LB simulations 0,0- 

To further test the accuracy of our method, we mea- 
sured the Stokes drag Fd on the colloid for various val- 
ues of the sphere radius a cs and gravitational field g. By 
varying the box size we find excellent agreement with an- 
alytic finite-size corrections [l^, from which the infinite 
box-size limit extrapolates to Fd — 71* s — iirrjfavs, as ex- 
pected for slip boundary conditions For the largest 
box sizes we compared the full velocity field around a 
single colloid to the known analytic result 0, and varied 
the ratio ao/a, finding that errors scale roughly linearly 
with this parameter. We choose ao/a w |, which leads to 
a relative error in the full velocity field of about 2%, sim- 
ilar to what is used in LB 8, 9], and sufficiently accurate 
for the kinds of questions we investigate 

The sedimentation runs were performed in a periodic 
box of dimensions L x = L y = 32ao and L z — 96ao, with 
N = 8 to 800 colloids and N f » 5 x 10 5 SRD parti- 
cles, corresponding to an average of about 5 particles per 
coarse- graining cell volume a^. The system size is simi- 
lar to some successful LB simulations |8|. A gravitational 
field g, applied in the z direction, was varied to produce 
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FIG. 1: Average sedimentation velocity, v s normalized by the 
Stokes velocity v®, as a function of volume fraction <f> for vari- 
ous Peclet numbers. The reduction of v s due to hydrodynamic 
backflow effects is independent of Pe. Dashed lines corre- 
spond to two versions of the semi-empirical Richardson-Zaki 
The dotted line is another theoretical 
prediction taking higher order HI into account [l^|. Ignoring 
hydrodynamics leads to Vs/v® = 1 — (dash-dotted line). 

different Pe. The simulations were run from 200 Stokes 
times t s (for Pe = 0.1) to 30, 000 t s (for Pe = 15), where 
ts = a/v s is the time it takes a sphere to sediment one 
particle radius. We verified that there was no drift in 
averages after about 100 ts, so that the suspension is in 
steady state. To check that our system is large enough, 
we performed some runs for double the box size described 
above, as well as for ao/a = 4, finding no significant 
changes in our conclusions 

The average sedimentation velocity v s for different 
Peclet numbers and different sphere packing fractions 
<j> = ^irpa 3 , with p the colloid number density, is shown 
in Fig. n At low densities the results are consistent 
with the Batchelor lawQ, while at higher densities they 
compare well with a number of other forms also de- 
rived for the Pe — > 00 limit 0). Although one might 
naively expect that the effect of HI becomes weaker for 
Pe < 1, we observe that the results for all Peclet numbers 
0.1 < Pe < 15 lie exactly on the same curve. Taking into 
account only Brownian fluctuations gives v s = v®(l — <f>) 
P|, which heavily underestimates backflow effects. This 
is strong evidence that purely hydrodynamic arguments 
are still valid in an average sense at low Pe. 

We next discuss velocity fluctuations around the av- 
erage. In colloidal systems the instantaneous velocity 
fluctuations 5w = v — v s are dominated by thermal fluc- 
tuations, with a magnitude determined by equipartition: 
A?4 = k B T/m. Here m is the mass of a colloid. To 
disentangle the hydrodynamic fluctuations from thermal 
fluctuations, we describe spatial and temporal correla- 
tions in the velocity fluctuations. The spatial correlation 
of the z component (parallel to the sedimentation) of the 
velocity fluctuations can be defined as 

C z {v) = (Sv z (0)Sv z (r)} , (2) 
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FIG. 2: Spatial correlation functions of the parallel (z) com- 
ponent of the velocity fluctuations as a function of distance 
perpendicular (a) and parallel (b) to the external field, for 
three different volume fractions (0 = 0.02 (grey symbols), 
(j> = 0.04 (white), <f> = 0.086 (black)) and different Peclet 
numbers. The correlation functions are scaled with v 2 to 
emphasize hydrodynamic fluctuations. The insets show how 
C z (r), scaled with C(0) = fcsT/m, increases with Pe. 

where (...) represents a time average over many parti- 
cles. The distance vector r is taken perpendicular to sed- 
imentation, C z {x), or parallel to it, C z (z). Similarly, the 
temporal correlation of the z component of the velocity 
fluctuations can be defined as 



C z (t) = (8v z (0)Sv z (t)) 



(3) 



In Fig. [2] we plot C z (r) , which shows a positive spa- 
tial correlation along the direction of flow, and an anti- 
correlation perpendicular to the flow, very much like that 
observed in experiments^). The inset of Fig.^a) shows 
that at Pe=l the correlation in the perpendicular di- 
rection, C z (x), is almost negligible compared with the 
thermal fluctuation strength fcsT/m, whereas for larger 
Pe, distinct regions of negative amplitude emerge, which 
grow with increasing Pe. Similarly, the inset of Fig. Elb) 
shows correlations in the parallel direction that rapidly 
increase with Pe. For the highest Peclet numbers studied 
(5 < Pe < 15), the amplitudes of these correlations grow 
proportionally to v 2 , as shown in the main plots of Fig-El 
Unfortunately, because the division by amplifies the 
statistical noise, we are unable to verify whether this scal- 
ing persists for Pe < 5. The minimum in Fig.^Ja) is lim- 
ited by the box size. We checked this by simulating larger 
systems: the correlation size increased linearly with box 
dimensions 0, as found for suggesting that the 

hydrodynamic velocity fluctuations are unscreened. 

The following timescales are important for temporal 
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FIG. 3: Temporal correlation functions of the z component 
of the velocity fluctuations for <f> = 0.02 and different Peclet 
numbers, (a) Time is scaled with the Brownian relaxation 
time tb — m/j and the velocities are scaled with the thermal 
fluctuation strength fcsT/m. The straight line is the hydro- 
dynamic long time tail Bt~' i/2 with B' 1 = 12p f k B T(-Kvf /2 
|15| . The results for Pe < 1 are indistinguishable, (b) Time is 
scaled with the Stokes time ts = a/v s and the velocities are 
scaled with v 2 to highlight hydrodynamic velocity fluctua- 
tions. The straight line is a fit demonstrating the exponential 
decay of non-equilibrium hydrodynamic fluctuations. 

correlations: In a liquid, the solvent relaxation time 
Tf, typically of order 10 _14 s[l|, is the smallest relevant 
timescale. A (larger) Brownian particle experiences ran- 
dom forces and a friction 7. As a consequence, it loses 
memory of its initial velocity after a time tb ~ rn/^f, 
which is typically of the order of 1CP 9 s 1]. For time 
scales larger than tb, the particle experiences diffusive 
behavior and after a time T£> = a /D 3> tb it has trav- 
eled over its own radius. For correct coarse-grained tem- 
poral behavior, the timescales don't need to be identical 
to the underlying fluid, but it is important that they are 
clearly separated 9]. This is indeed the case for our choice 
of SRD parameters, where r/ « At c = 0.1, tb ~ 2, and 
td ~ 200 (in units ao{mf /ksT)^). Since the Stokes time 
ts = a/v s = rn/Pe must be ^> tb, this this sets a limit 
on the maximum Pe number for these parameters. 

Fig. |31 shows the temporal correlation functions along 
the direction of sedimentation. At short times the 
behavior is well described by exponential Brownian 
relaxation^: C s hort(t)=Av^ exp (— £/tb). At intermedi- 
ate times it follows the well known algebraic long time tail 
Ciong(t)=Bt~ 3 / 2 , associated with the fact that momen- 
tum fluctuations diffuse away at a finite rate determined 
by the kinematic viscosity v. Analytical equilibrium cal- 
culations of .B|l5j exactly fit the low Pe nonequilibrium 
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FIG. 4: Scaling of the hydrodynamic relaxation times (left 
scale) and velocity fluctuation amplitudes (right scale) with 
volume fraction. Straight lines are expected scalings for an 
unscreened svstem|4l ll6|l . 

results in Fig. |3{a) with no adjustable parameters! 

Several experimental studies[l4| on the sedimentation 
of non-Brownian (Pe — > go) particles have found an ex- 
ponential relaxation of the form 

C z (t) = Av 2 H exv{-t/T H ). (4) 

This non-equilibrium hydrodynamic effect takes place 
over much longer time-scales than the initial exponential 
Brownian relaxation. The double-logarithmic Fig. EJa) 
shows that a new mode of fluctuations becomes distin- 
guishable in our simulations for Pe > 1. In Fig- Etb) the 
correlation functions are scaled with v 2 s to highlight the 
nonequilibrium hydrodynamic fluctuations. For Pe > 10 
the fluctuations scale onto a single exponential master 
curve, similar to the high-Pe experiments whereas 
for lower Pe deviations are seen. From the exponential 
fit to Eq. we can estimate the relaxation time th and 
the amplitude Av 2 H of the hydrodynamic fluctuations. 
These are shown in Fig.0|for different volume fractions <fi. 
The scalings of the relaxation time and fluctuation am- 
plitude with (p are consistent with Av 2 H « v 2 L<p/ (^na) 
and rfj « L 2 / Av 2 H « tg^irL/ '(4>a), predicted for un- 
screened hydrodynamic fluctuations by a simple heuristic 
argument akin to that used by Caflish and Luke Q. 

As seen in Fig. [3] the short time velocity fluctuations 
are dominated by thermal fluctuations at all Peclet num- 
bers studied. By comparing Avh with Avt, we estimate 
the critical Pe*, above which hydrodynamic fluctuations 
are larger than thermal fluctuations for all t [l7| : 

Av H (k B TL<j>p c )* ^ _ Pe 

For example, for polystyrene colloids in water (77 = 
1CT 3 Pa s, T = 300 K, p c = 1050 kg m~ 3 ), Pe* « 

[(a/10- 14 m)/((j)L/a)] K For = 0.001, a = 10~ 6 m, and 
L/a = 100 (smaller than the screening length at this con- 
centration 0,0), we find a large value: Pe* «3x 10 4 . 

In conclusion, we have adapted a mesoscopic simula- 
tion method, SRD0, to study sedimentation at finite 
Peclet numbers. Hydrodynamic backflow corrections re- 
duce the average sedimentation velocity v s , irrespective 



of Pe. Thus, even when HI are relatively small, Brownian 
dynamics simulation methodsQ] will yield qualitatively 
incorrect results for this problem. Long-time nonequi- 
librium velocity fluctuations become evident for Pe > 1, 
and scale like those for Pe — > 00, while short time fluctu- 
ations are dominated by Brownian forces up to surpris- 
ingly large Pe. In other words, neither hydrodynamic 
interactions nor Brownian forces can be ignored for a 
significant parameter regime. 
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